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Abstract 



We propose a method for solving boundary value and eigenvalue problems for the 
elliptic operator D — divp grad +q in the plane using pseudoanalytic function theory and 
t-H in particular pseudoanalytic formal powers. Under certain conditions on the coefficients 

P and q with the aid of pseudoanalytic function theory a complete system of null solutions 
of the operator can be constructed following a simple algorithm consisting in recursive 
integration. This system of solutions is used for solving boundary value and spectral 
problems for the operator D in bounded simply connected domains. We study theoretical 
and numerical aspects of the method. 

X 

1 Introduction 

The main numerical techniques for solving problems related to elliptic linear partial differ- 
ential equations with variable coefficients in one way or another involve a discretization of a 
domain and solution of systems of thousands of algebraic equations. Seldom the method of 
separation of variables is applied due to its natural limitations related to the requirements 
of a complete agreement between the geometry of the domain and the symmetry of the coef- 
ficients. Moreover, the method of separation of variables implies solution of Sturm-Liouville 
spectral problems which is not an easy task itself. 

In the present paper we propose a different method based on some old and new results 
from pseudoanalytic function theory [3], |15j . Its applicability is not so universal as the 
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wishes to thank support from the CONACYT and the National Polytechnic Institute for the possibility of a 
postdoctoral stay in the Department of Mathematics of the CINVESTAV in Queretaro, as well as from the 
SIBE program of the National Polytechnic Institute, Mexico. 
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applicability of the finite difference method or the finite element method. First of all, it is 
applicable to problems in bounded domains in the plane and up to now only for the operator 
divpgrad+g. Moreover, at present we can apply the method only when the equation 



nonvanishing in the domain of interest and representable in the form / = S(s)T(t) where 
s and t represent an orthogonal coordinate system. This separable form of / may cause 
associations with the method of separation of variables. Nevertheless this is a completely 
different technique, based on different ideas and free of the mentioned above limitations of 
the method of separation of variables. 

The heart of the method is the construction of a complete system of solutions for ([!]) in 
the domain of interest, complete in the sense explained below (see Sect. 1.3] for related 
ideas and additional details). The system of solutions is used for approximating the solution 
of a boundary value problem. Due to the linearity of equation ([!]) after the construction of 
the system of solutions the problem reduces to approximation of boundary conditions, for 
which a variety of methods can be used. Here we apply the collocation method. 

The complete system of solutions is constructed in the following way. The knowledge 
of a particular solution of ([I]) allows us to propose a corresponding Vekua equation [19] 
closely related to |IJ in the sense that the real part of any of its solutions has the form p l l 2 u 
where u is a solution of ([I]), and vice versa given u one can easily construct a corresponding 
solution of the Vekua equation [12] , [15] . The relation between ([l]) and the Vekua equation 
is similar to the relation between the Laplace equation and the Cauchy-Riemann system. 
L. Bers developed [3J, [I] a theory of so-called pseudoanalytic formal powers. They are 
generalizations of the analytic powers (z — zq) 11 in the sense that they are solutions of the 
corresponding Vekua equation and behave locally like the analytic powers. The theory of Bers 
includes generalizations of Taylor series, Runge's theorem and other basic facts from analytic 
function theory. Thus, under certain quite natural conditions the system of pseudoanalytic 
formal powers is complete in the space of all pseudoanalytic functions (solutions of the 
Vekua equation) in the same sense as the system of powers (z — zo) n is complete in the 
space of analytic functions. To construct the pseudoanalytic formal powers the knowledge 
of a corresponding generating sequence is required. Recently [H], [15] an algorithm for 
construction of generating sequences under additional conditions on the coefficients in the 
Vekua equation was proposed. This implies that when / = p^^uq is representable in a 
separable form the complete system of formal powers for the Vekua equation associated with 
([!]) can be constructed explicitly following Bers' recursive procedure. 

We investigate the efficiency of the proposed method which we call MPFP, the Method 
of Pseudoanalytic Formal Powers. We show its fast convergence and compare its accuracy 
with that of the finite element method. In general, we show that in problems addmitting 
the explicit construction of formal powers and hence the application of the MPFP its use 
is advantageous compared to other computational techniques based on discretization of the 
problem. 

It is worth mentioning that the MPFP is a direct generalization of the method of har- 
monic polynomials for solving boundary value problems for the Laplace equation which 
has been considered in dozens of works (see, e.g., [5], [TO], [H], [IB]). Indeed, in a spe- 
cial case when p = 1, q = and Uq = 1 the corresponding complete system of solutions 
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constructed by means of the MPFP coincides with the system of harmonic polynomials 

{Re(z-z ) n , Im(z - z ) n }Zo- 

The knowledge of a complete system of solutions for an equation corresponding to any 
value of a spectral parameter allows one to use it for solving eigenvalue problems. We consider 
this possibility in section|6| The numerical results are highly promising, and it is clear that in 
the case of eigenvalue problems as well as for boundary value problems further work should 
be done in investigation of optimal ways of application of the MPFP. For example, for solving 
eigenvalue problems by means of the MPFP we used the simplest possible idea reducing the 
problem to calculation of zeros of a certain determinant obtained by evaluating the first N 
solutions from the constructed complete system in N points on the boundary of the domain 
under consideration. Meanwhile, in principle, this natural approach works there exist other 
techniques offering different ways of using the available exact solution systems (see [21 Sect. 
1.13], where similar questions are discussed). 



2 Factorization of the operator divj9grad+g. 

Let 0, be a domain in R 2 . Throughout the whole paper we suppose that £1 is a simply 
connected domain. Denote dz = \ + and d z = \ ( — iJ0 . By C we denote the 
operator of complex conjugation. 

Note that the operator dz applied to a real valued function ip can be regarded as a kind 
of gradient, and if we know that dz<p = $ in a whole complex plane or in a convex domain, 
where = 3>i + i&2 is a given complex valued function such that its real part <J?i and 
imaginary part $2 satisfy the equation 

dy$X ~ d x $ 2 = 0, (2) 

then we can reconstruct (p up to an arbitrary real constant c in the following way 

ip(x,y) = 2(f $ 1 (ri,y)dri+ f $ 2 (xo,O d U + c 
\Jx J yo / 

where (xo>2/o) is an arbitrary fixed point in the domain of interest. Note that this formula 
can be easily extended to any simply connected domain by considering the integral along an 
arbitrary rectifiable curve T leading from (xo,yo) to (x,y) 

<p(x,y) = 2(J^$ 1 dx + $ 2 dy\+c. (3) 

By A we denote the integral operator in (|3|: 

A[§](x,y) = 2 ( f <$> l {ri,y)dr,+ f $ 2 (x , c. 

\Jx Jyo J 

Thus if $ satisfies Q, there exists a family of real valued functions cp such that d&p = 
given by the formula <p = A[$]. 

The following result is in the core of the method proposed in the present work. 
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Theorem 1 Let p and q be real valued functions, p G C 2 (Q) and p ^ in £1, uq be a 
positive particular solution of the equation 



(divp grad +q)u = in Q. (4) 

Then for any real valued continuously twice differentiable function ip the following equality 
holds 

- A {d^pgT^+q)^=p 1 ' 2 {d z + f -jC^ fdg-^Cjp 1 ' 2 ?, (5) 

where 

f=p 1/2 u - (6) 
Remark 2 Let q = 0. Then uq can be chosen as uq = 1. Hence |5[) gives us i/ie equality 

Let / be a real function of x and y. Consider the Vekua equation 

Wz = jW in n. (7) 

This equation plays a crucial role in all that follows and hence we will call it the main Vekua 
equation. We notice that the operator of this equation is precisely the second factor in ([5]). 
Denote W\ = ReW and W 2 = ImW. 

Theorem 3 [13] Let W = W\ + iW% be a solution of (71). Assume that f = p l / 2 UQ, where 
uq is a positive solution of in Q. Then u = p~ x I^V\ is a solution of in 0,, and 
v = p x l 2 W2 is a solution of the equation 

(div - grad +q\)v = in O, (8) 
P 



where 



p \p \ p Uq \ Uq 



Theorem [3] shows us that as much as real and imaginary parts of a complex analytic 
function are harmonic functions, the real and imaginary parts of a solution of the main 
Vekua equation |7| multiplied by p~ 1 / 2 and p 1 / 2 respectively are solutions of the associated 
elliptic equations Q and The following natural question arises then. We know that 
given an arbitrary real valued harmonic function in a simply connected domain, a conjugate 
harmonic function can be constructed explicitly such that the obtained couple of harmonic 
functions represent the real and imaginary parts of a complex analytic function. What is 
the corresponding more general fact for solutions of associated elliptic equations Q and ^ 
(which we slightly generalizing the definition of I. N. Vekua call metaharmonic functions). 
The precise result is given in the following theorem. 
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Theorem 4 [13] Let f = p 1 ^ 2 ^, where uq is a positive solution of in a simply connected 
domain Q and u be a solution of Q). Then a solution v of with q\ defined by |5p such 
that W = p l l 2 u + ip~ x l 2 v is a solution of ^), is constructed according to the formula 

v = u 1 A(ipulc\{u G l u)). (10) 

Let v be a solution of Q), then the corresponding solution u of such that W = p l / 2 u + 



vp 



-1/2 



V IS 



a solution of ([71), is constructed according to the formula 
u = -u A(ip~ 1 UQ 2 c\(u v)). 



(11) 



Remark 5 When p = 1, q = and uq = 1, equalities (10) and (11) turn into the well 
known formulas in complex analysis for constructing conjugate harmonic functions. 



3 Formal powers 

Briefly speaking formal powers are solutions of a Vekua equation 

Wz = aW + bW (12) 

(with a and b being complex valued functions) generalizing the usual analytic powers {(z — zo) n } 
in the sense that locally when z — > zq they behave asymptotically like the usual powers and 
under some additional conditions on the coefficients a and b they form a complete system 
in the space of all solutions of the Vekua equation in the same sense as the analytic powers 
{(z — zo) ri }^L l° rm a complete system in the space of analytic functions. Generalizations 
of the extension theorem, the Runge theorem and of other important results about the con- 
vergence of corresponding series are valid. The construction of formal powers is one of the 
main problems of pseudoanalytic function theory. Recently it was solved [13] , |15j for a wide 
class of Vekua equations of the form ([7]) which as was shown in the preceding section are of 
main interest for studying problems for second order equations of the form Q. 

The main ingredient for obtaining the explicit form of formal powers for a certain Vekua 
equation is the generating sequence, a concept introduced by Bers. If one knows a generating 
sequence for a given Vekua equation then the construction of formal powers reduces to a 
simple algorithm. Here we briefly explain the main ideas and steps refering the reader to [3 j 
and |15j for further details. 



3.1 Generating pair and generating sequence 



Definition 6 A pair of solutions F and G of a Vekua equation (12) in £1 possessing partial 
derivatives with respect to the real variables x and y is said to be a generating pair if it 
satisfies the inequality 

Im(FG) > inn. (13) 



Condition (13) implies that every complex function W defined in a subdomain of 17 
admits the unique representation W = 4>F + ipG where the functions (ft an d "0 are rea l 
valued. Thus, the pair (F,G) generalizes the pair which corresponds to usual complex 
analytic function theory. The following expressions are known as characteristic coefficients 
of the pair (F, G) 
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FGz - FzG _ FGz - FzG 

Hm - - fc-FG_ ' b(F > G) ~ FG-FG ' 

. FG Z - F Z G FG 7 - F Z G 

Aiprr, = = — = — • -D 



™ " " pa-PC ' °^ G) ~ FG-FG 



If (F,G) is a generating pair of a Vekua equation (12) then clif,g) = a an d brp,G) = b. 
The other two characteristic coefficients are related to the concept of a derivative [3j . The 

di F ,G)W 
dz 

form 



(F, G)-derivative W = (F ff — of a continuously differentiable function W exists and has the 



W = W z -A {F>G) W-B (FjG) W (14) 

if and only if 

Wz = a(F,G)W + \f,g) w - 
Solutions of this equation are called (F, G)-pseudoanalytic functions. 

Definition 7 Let (F,G) and (Fi,Gi) - be two generating pairs in 0. (Fi,G\) is called 
successor of (F, G) and (F, G) is called predecessor of (Fi, G±) if 

a (F\,G\) = a (F,G) and \F\,Gx) = —B(F,G)- ( 15 ) 

This definition arises naturally in relation to the notion of the (F, G)-derivative due to 
the following fact. 

Theorem 8 Let W be an (F, G)-pseudoanalytic function and let (Fi,G\) be a successor of 
(F,G). Then W is an (F\,G{)-pseudoanalytic function. 

Thus, to the difference of analytic functions whose derivatives are again analytic, the 
(F, G)-derivatives of pseudoanalytic functions are in general solutions of another Vekua equa- 
tion with the coefficients given by (15). Obviously this process of construction of new Vekua 
equations associated with the previous ones via relations ( 15 ) can be continued and we arrive 
at the following definition. 

Definition 9 A sequence of generating pairs {(F m ,G m )} , m = 0,±1,±2, . . is called a 
generating sequence if (F m+ i, G m+ \) is a successor of (F m ,G m ). If (Fq,Go) = (F,G), we 
say that (F,G) is embedded in {(F m ,G m )}. 

Definition 10 A generating sequence {(F m , G m )} is said to have period \i > if(F m+ ^, G m+At ) 
is equivalent to (F m ,G m ) that is their characteristic coefficients coincide. 

We will need the following notation introduced by Bers. The (F, G)-integral is defined 
as follows 

f f 2£r f 2F 

where T is a rectifiable curve leading from zq to z%. 
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Let W be an (F, G)-pseudoanalytic function. Using a generating sequence in which (F, G) 
is embedded we can define the higher derivatives of W by the recursion formula 

W M =W . W [m + l ] = d iF m ,G )W m = Q1 

dz 

A generating sequence defines an infinite sequence of Vekua equations. If for a given 
(original) Vekua equation we know not only a corresponding generating pair but the whole 
generating sequence, that is a couple of exact and independent solutions for each of the 
Vekua equations from the infinite sequence of equations corresponding to the original one, 
we are able to construct an infinite system of solutions of the original Vekua equation as is 
shown in the next definition. 

Definition 11 The formal power Zm\a,zo;z) with center at zq 6 f2, coefficient a and 
exponent is defined as the linear combination of the generators F m , G m with real constant 
coefficients A, \i chosen so that \F m {zQ) + fiG m (zo) = a. The formal powers with exponents 
n = 1,2,... are defined by the recursion formula 

Z$(a,z Q \z) = n I Z%~t\a,z ;()d(F m ,G m )C- (16) 

J z 

This definition implies the following properties. 

1. z£\a , zo; z) is an (F m , G m )-pseudoanalytic function of z. 

2. If a' and a" are real constants, then Zm\a'+ia", zq; z) = a'Zm 1 (1, zq; z)+a" Z%, (i, zq; z). 

3. The formal powers satisfy the differential relations 

d(F m ,G m )Z^(a,z ;z) _ ( w -i), , 

JZ - nZ m+l [ a i z 0,Z). 

4. The asymptotic formulas 

Zj$(a,zo;z) ~ a(z- z ) n , z -> z (17) 

hold. 

Assume now that 

oo 

W(z) = Y,Z (n) (a n ,z ;z) (18) 

n=0 

where the absence of the subindex m means that all the formal powers correspond to the same 
generating pair (F,G), and the series converges uniformly in some neighborhood of zq. It 
can be shown that the uniform limit of pseudoanalytic functions is pseudoanalytic, and that 
a uniformly convergent series of (F, G)-pseudoanalytic functions can be (F, G)-differentiated 
term by term. Hence the function W in (18) is (F, G)-pseudoanalytic and its rth derivative 
admits the expansion 

oo 

Wtt(z) = £n(n - 1) • • • (n - r + l)Z^~ r \a n , z ; z). 
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From this the Taylor formulas for the coefficients are obtained 

a. = ^>. (19) 
n! 

Definition 12 Let W(z) be a given (F, G)-pseudoanalytic function defined for small values 
of \z — zq\. The series 



J2z^Ha n ,z ;z) (20) 



n=0 



with the coefficients given by (19) is called the Taylor series ofW at zq, formed with formal 
powers. 

The Taylor series always represents the function asymptotically: 

N 

W(z)-^2z^(a n ,z ;z) = o{\z-z \ N+1 ) , z -> z , (21) 

n=0 

for all N. This implies (since a pseudoanalytic function can not have a zero of arbitrarily high 
order without vanishing identically) that the sequence of derivatives {VF^ (zo)} determines 
the function W uniquely. 



If the series (20) converges uniformly in a neighborhood of zo, it converges to the function 

W. 

3.2 Convergence theorems 

S. Agmon and L. Bers pQ and L. Bers developed a theory of expansions in pseudoanalytic 
formal powers which in its generality is presented in [3], [3]. We do need here the general 



results concerning a general Vekua equation (12). Fortunately the situation with the main 



Vekua equation ([7]) in a bounded simply connected domain under quite natural conditions on 
the function / is much easier than in the general case, and we have the following expansion 
theorem and Runge theorem |14j . |15j . 



Theorem 13 Let D be a disk of a finite radius R and center zq, and f 6 C 1 (-D) be positive 
in D. Then any solution W of ([?[) in D admits a unique normally convergent expansioi^ 
of the form W(z) = ZZo Z {n) (a n , z ; z). 

Theorem 14 any solution W of ([?p defined in a simply connected domain can be expanded 
into a normally convergent series of formal polynomials (linear combinations of formal pow- 
ers with positive exponents). 

Remark 15 This theorem admits a direct generalization onto the case of a multiply con- 
nected domain (see 

We mention here another important result obtained by Menke in |17j which gives a useful 
estimate for the rate of convergence of the series from the preceding theorem in the case when 
W is a Holder continuous function up to the boundary of the domain of interest. 



1 Following [3], [5] we shall say that a sequence of functions W n converges normally in a domain f2 if it 
converges uniformly on every bounded closed subdomain of Q. 
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Theorem 16 Let W be a pseudoanalytic function in a domain f2 bounded by a Jordan curve 
and satisfy the Holder condition on dfl with the exponent a (0 < a < 1). Then for any e > 
and any natural n there exists a pseudopolynomial of order n satisfying the inequality 

\W(z)-P n (z)\< C ^ for any z eTl 
n a e 

where the constant does not depend on n, but only on e. 

The following statements are direct corollaries of the relations established in section [2] 
between pseudoanalytic functions (solutions of 0) and solutions of second-order elliptic 
equations, and of the convergence theorems formulated above. Here we assume the existence 
of a positive solution uq of Q^in the domain Q and the function / in Q to be defined by 
/ = p 1 I 2 uq and belong to C'^Jtl). 

Definition 17 Let u(z) be a given solution of the equation defined for small values of 
\z — zq\, and let W(z) be a solution of |?p constructed according to theorem^ such that 
ReW = p l / 2 u. The series 

oo 

p- 1 / 2 (z)^2ReZ^(a n ,z ;z) (22) 

n=0 

with the coefficients given by (jg[ ) is called the Taylor series of u at zq, formed with formal 
powers. 

Theorem 18 f!3\/ . Ji5]/ Let u{z) be a solution of defined for \ z — zq\ < R. Then it admits 
a unique expansion of the form 



u(z)=p- l / 2 (z)J2^Z {n) (a n ,z ;z) 

n=0 

which converges normally for \z — zq\ < R. 

Theorem 19 An arbitrary solution of defined in a simply connected domain where there 
exists a positive particular solution u$ such that f = p l l 2 u$ 6 C 1 (il) can be expanded into a 
normally convergent series of formal polynomials multiplied byp~ l l 2 . 

More precisely the last theorem has the following meaning. Due to Property 2 of for- 
mal powers we have that Z^ n \a, zq; z) for any Taylor coefficient a can be expressed through 



Z( n \l, zq; z) and Z^ n \i, zq; z). Then due to theorem 14 any solution W of (j7|) can be 
expanded into a normally convergent series of linear combinations of Z^ n \l, zq; z) and 
Z^ n \i, zq; z). Consequently, any solution of Q can be expanded into a normally conver- 
gent series of linear combinations of real parts of Z^ n \l, zq; z) and Z^ (i, zq; z) multiplied 
by p~ x l 2 . 

Obviously, for solutions of Q the results on the interpolation and on the degree of 
approximation like, e.g., theorem 16 are also valid. 

Let us stress that theorem 19 gives us the following result. The functions 

\p-^ 2 (z) Re ZW(1, z ; z), p- l ' 2 (z) Re zW(i, z ; z)}°° (23) 

^ J 71=0 
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represent a complete system of solutions of Q in the sense that any solution of Q can be 
represented by a normally convergent series formed by functions ( 23 ) in any simply connected 
domain fi where a positive solution of Q exists, and the rate of convergence of the series 
can be estimated with the aid of theorem f" 



3.3 Explicit construction of generating sequences and formal powers 

The results of section [2] show us that the theory of the elliptic equation 

(divpgrad +q)u = 

is closely related to equation ([7]): 

W, = jW. (24) 
It is interesting that for this equation we always know a generating pair. Namely, it is easy 



to see that the functions F = f and G = j satisfy (24) together with the condition (13) 
Then the corresponding characteristic coefficients ^4(_f,G) an d Brp t Q\ have the form 

fz 

A (F,G) = 0) B (F,G) = y 5 



and the (F, G)-derivative according to ( 14 ) is defined as follows 



W = W r - hw = \d 7 - he \ W. 



f \ f 
Due to Theorem [8] we obtain the following statement 



Proposition 20 Let W be a solution of (24). Then its (F,G)- derivative, the function w 
W is a solution of the equation (c\+ ^fC\ w = 0. 



In spite of having given a generating pair for (24) in general it is not known how to 
construct a corresponding generating sequence necessary for calculating the system of formal 
powers. Nevertheless a recent result from |14j . |15j which we formulate in the following 
statement gives an answer to this question in a quite general situation. 

Theorem 21 Let F = S(s)T(t) and G = g^j'^ where S and T are arbitrary differ entiable 
nonvanishing real valued functions, $ = s+it is an analytic function of the variable z = x+iy 
in such that & z is bounded and has no zeros in f2. Then the generating pair (F,G) is 
embedded in the generating sequence (F m , G m ), m = 0, ±1, ±2, ... in £1 defined as follows 

F m = (& z ) m F and G m = (& z ) m G for even m 

and 

F m = ^— ^ — F and G m = ($ 2 ) m S 2 G for odd m. 
S z 
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In order to appreciate the generality of this construction let us remind that orthogonal 
coordinate systems in a plane are obtained (see |16| ) from Cartesian coordinates x, y by 
means of the relation 

s + it = <&(x + iy) 

where $ is an arbitrary analytic function. Quite often a transition to more general coordi- 
nates is useful 

£ = £( s ), V = V (t). 

£ and r/ preserve the property of orthogonality. To illustrate the point, besides the obvious 
example of Cartesian coordinates which are generated by the analytic function z we give 
some other examples taken from [16J. 

Example 22 Polar coordinates 

s + it = ln(x + iy), 

s = In \/ x 2 + y 2 , t = arctan — . (25) 

x 

Usually the following new coordinates are introduced 

r = e s = \/ x 2 + y 2 , <p = t = arctan — . 

x 



Example 23 Parabolic coordinates 

s + it 



y/2 

s = \/r + x, t = y/r — x. 
More frequently the parabolic coordinates are introduced as follows 

d = s 2 , n = t 2 . 

Example 24 Elliptic coordinates 

. x + iy 

s + it = arcsm , 

a 

s\ - s 2 . , Sl + s 2 

sm s = , cosn t = 

2a 2a 

where s± = yjx + a) 2 + y 2 , s 2 = yjx — a) 2 + y 2 . The substitution 

£ = sins, T) = cosht 

is frequently used. 
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Example 25 Bipolar coordinates 



tanh s 



s + it = In 



2ax 



a + x + iy 
a — x — iy' 



a 2 + x 2 + y 2 ' 
T/ie following substitution is frequently used 



tant 



lay 



a 2 — x 2 — y 2 



= e S , rj = ir — t. 

The last theorem opens the way for explicit construction of formal powers corresponding 



to the main Vekua equation ( 24 ) in the case when / has the form 

/ = S(s)T(t) 



(26) 



and hence for explicit construction of complete systems of solutions for corresponding second- 
order elliptic equations admitting a particular solution of this form. 



4 Description of the method 

We consider boundary value problems of Dirichlet, Neumann or mixed type for the elliptic 
equation of the form Q in a bounded, simply connected domain Qct 2 . The main assump- 
tion required for the applicability of the method of pseudoanalytic formal powers (MPFP) 
is the existence in Q of a positive solution uq such that the function / = p 1 / 2 uo £ C 1 (J7) 



be representable in a separable form (26) in an orthogonal coordinate system. Let us stress 



that very often such a particular solution uq is readily available. The simplest example 



of such situation is when q = and p is of the form (26). For example, the cases when 
p(x,y) = X(x)Y(y) or p = p(^x 2 + y 2 ) frequently occur in practice [8]. 
When the equation of the form 

(-A + q{y))u{x,y) + X 2 u(x,y) =0 (27) 

is considered, it is sufficient to obtain a particular solution for the ordinary differential 
equation 

(--^2+q(y)My) = o. (28) 

Then a particular solution of ( |27| ) can be constructed as follows 

u (x,y) = e Xx h(y). (29) 

It has a convenient separable form. Notice that in this example we come to an important 
open problem. It is related to the requirement that u$ should be different from zero in 
the domain of interest. Meanwhile in many practically significant situations it is easy to 
guarantee that h(y) ^ when (x, y) £ O, sometimes this condition becomes a considerable 



obstacle. Moreover, when A in (27) is purely imaginary, the solution (29) is not acceptable 



because it is no longer real valued. In this case one should take instead of e lkx , where A = ik 
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a solution in the form of sin kx or cos kx but then for a big domain or large k one cannot 
avoid the appearance of zeros of the resulting particular solution no and in this case the 
proposed scheme in general does not work. 

One possibility to overcome this problem is to include under consideration a complex 
valued particular solution no but then we would need to consider a corresponding bicomplex 
main Vekua equation (see |13| and [E]). Then the whole algorithm for the construction of 
generating sequences and formal powers would go through with no modification compared 
to the complex case, but up to now there is no proof of the completeness of the system of 
formal powers for a bicomplex Vekua equation. As a consequence there is no guarantee that 



the infinite system of exact solutions obtained similar to (23) will be complete in the space 
of solutions of Q in Q. Our conjecture is that at least in the case when Uq(x, y) = g(x)h(y) 
where g and h are complex valued nonvanishing functions the system of formal powers for 
the corresponding bicomplex main Vekua equation is complete in the same sense as was 
established earlier for the complex case. We continue this discussion in section [6] where we 



use complex valued particular solutions of the form ( 29 ) for solving eigenvalue problems for 
operators of the form —A + q(y). 

Turning back to equation Q we assume that it admits a positive solution no in the 
domain f2 such that / = p 1//2 no 6 C 1 (0) is representable in a separable form (26) and that 



$ = s + it is an analytic function of the variable z = x + iy in Q such that $ 2 is bounded and 
has no zeros in Q. Then applying theorem |21| one can construct a corresponding generating 
sequence. Construction of formal powers {ZW(l,zo-,z), Z( n \i,z ;z)}™ =0 reduces then to 
the recursive algorithm described in Definition 1 1 , and in this way one obtains the complete 
system of solutions for Q in f2 given by |23| ). By construction Re Z^°\i, zq; z) = 0. Taking 
this into account we introduce the notations 

u 1 (z)=p- 1 / 2 (z)ReZW(l,z ;z), u 2 (z) = p-^\z) ReZ^(i, z ; z), 
u 3 (z)=p- 1 / 2 (z)ReZ ( - 2 \l,z ;z), u A (z) = p- 1/2 (z) ReZ^ 2 \i, z ;z), . . . 

and obtain the complete system of solutions for Q given by {no, ni, U2, ■ ■ ■}■ We look for an 
approximate solution of a boundary value problem for Q in the form 



N 

.N 



U 

k=0 



huk (30) 



where are real coefficients which should be found from boundary conditions. To obtain 
iV + i equations for finding {bk}^ =0 one can 
points Cj G dQ, we obtain N + 1 equations 



iV + 1 equations for finding {bk}^ =0 one can use, e.g., the collocation method. Chosing N + 1 



N 



Y / hB[u k ]{C j ) = v(C j ), j = 0,N 



k=0 

where v is a given function and B is the linear operator of the boundary condition. For 
the Dirichlet condition one has B[u] = u and for the Neumann condition, B[u] = - the 
normal derivative of u. Finding {bk}^ =0 we have the approximate solution u N . 

Thus, the proposed here MPFP belongs to the class of boundary methods because due 
to the linearity of the problem the function (30) is an exact solution of Q in and only 
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boundary conditions should be approximated. An estimate for the rate of convergence of 
the method is given in theorem 16 In the next section we discuss the numerical realization 
of MPFP and results of numerical tests. 



5 Approximate solution of boundary value problems 

As a first example we considered the Dirichlet problem for the equation 

(-A + c 2 )n = (31) 

where c is a real constant. The interest in this relatively simple equation is not due to its 
numerical simplicity. In fact this is not the case,- numerical solution of this equation is not 
less difficult than that of an equation with c being a reasonably good function with a range 
of values comparable with c. The attractiveness of this example consists in the possibility 
to calculate a large number of the functions Uk (see the preceding section) symbolically, 
using an appropriate software for symbolic calculations like Mathematica (Wolfram), Maple 
or Matlab. In this work we used Matlab 2006 and a PC of 2 GB in RAM and a processor 
of 1.73 GHz. Implementation of the symbolically calculated base functions Uk gives us the 
possibility to estimate the accuracy of the MPFP itself without considering the precision of 
recursive numerical integrations. We also compare the results obtained using symbolically 
calculated Uk with the results obtained purely numerically. 



For equation (31) it is easy to propose a positive particular solution. It can be chosen, 
e.g., as / = e cy . Then the first functions constructed as described in the preceding section 
taking as a center of the formal powers the origin will have the form [15J 



/ \ cy t \ cy i \ sinh(cy) 
uo{x, y) = e y , ui(x, y) = xe \ u 2 {x, y) = , 

(32) 

/ 2 y\ sinh(cy) 2xsinh(cy) 

U3{x,y) = [x z e cy H 2 — ' udx,y) = , •••• 

V c/ c c 

It is interesting to mention that using Matlab we obtained the first 101 functions of this sys- 



tem calculated symbolically. According to theorem 19 this system of solutions is complete in 
any bounded simply connected domain containing the origin. First we show results obtained 
with the help of the system of functions Uk calculated symbolically. 

5.1 Numerical results obtained with symbolically calculated base func- 
tions 

We begin with the unitary disk D with center in the origin. As a test exact solution we take 
the function 

u = e cx . (33) 



cx 



Thus, the problem we consider is to solve (31 ) in D with the boundary condition u\ an = e 
We look for an approximate solution u N in the form ( 30 ) with the base functions ( 32 ) . We 



use the collocation method for satisfying the boundary condition, the collocation points are 
distributed uniformly on dD. Their number is equal to the number of solutions Uk- 



14 



According to theory from subsection 3.2 the coefficients bj~ in (30) are obtained in the 



case under consideration from the Taylor coefficients which appear in (22). More precisely 



we have that according to theorem [18^ the solution u can be represented as follows 



ulz 



J2ReZ^ n \a n ,0;z) = ]T (a' n Re Z (n) (l, 0; z) + < Re Z^ n \i, 0; . 



n=0 



n=0 



where a n = a' n +ia'^ are the Taylor coefficients given by (19). In the case of the exact solution 



(33) the Taylor coefficients have the form [151 Sect. 7.3] 



°n = -r(l + i)- 



Thus, the exact values for the coefficients bk from (30) in our example are as follows 

c 2 

b Q = 1, b 1 = b 2 = c, 63 = 64 = —,.... 



Having compared the numerically calculated constants bk which we denote by bk for 
N = 34 with their exact values in the case c = 1 we obtained their coincidence up to 10~ 14 
for every k = 0, . . . , 34. For smaller values of c the situation is the sa me. The difference 
between bk and bk tends to become larger for larger values of c. In Table 5.1 we show results 
for c = 5 and N = 34. 

Table 1. Comparison of the values of bk and bk as k increases 



k 


The values of bk 


The values of bk 


bk ~ b k 


5 


20.83333333333382 


20.83333333333333 


0.00000000000048 


8 


26.04166666666448 


26.04166666666667 


0.00000000000219 


13 


15.50099206509864 


15.50099206349206 


0.00000000160657 


17 


5.38228885848900 


5.38228891093474 


0.00000005244574 


25 


0.19601580023149 


0.19603324996120 


0.00001744972971 


31 


0.00743227875004 


0.00729290364439 


0.00013937510565 


34 


0.00172611010091 


0.00214497166011 


0.00041886155920 



The maximum number of functions Uk that we used here is limited not by the possibility 
of obtaining them symbolically but rather by the time required for numerical calculations 
involving the corresponding quite long symbolic expressions. 

In the following two tables, the convergence of MPFP is shown by comparison of the 
maximum absolute error obtained for different values of N, for the case c = 1 and c = 5. 

Tables 2 and 3 Maximum absolute error depending on N for c = 1 and c = 5 



N 


Maximum absolute error 




N 


Maximum absolute error 


8 


0.00698626935341 




6 


3.59578971016677 x W 2 


14 


2.534633673767495 x 10" 5 




14 


22.38029523897584 


22 


1.432881036045330 x 10" 9 




22 


0.73431266884919 


28 


4.276579090856103 x lO" 13 




32 


0.00194275813006 


32 


1.776356839400251 x 10~ 15 




44 


0.59167057031573 x 10~ 7 


36 


8.881784197001252 x 10" 16 




54 


0.72795103278622 x lO" 11 


38 


1.110223024625157 x 10~ 15 




60 


8.781864124784988 x 10~ 14 
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In Table 4 for a fixed number (N = 34) of the base functions Uk we show the dependence 
of the maximum absolute error on the parameter c. Here we also indicate the maximum 
absolute error obtained for the same problem using the standard PDE tool of Matlab. 

Table 4. Performance of MPFP compared to Matlab's PDE tool in terms of the maximum 
absolute error for increasing values of c as N = 34 



c 


Maximum absolute error of MPFP 


PDE tool (2129 nodes) 


0.1 


0.89 x lO" 15 


1.5 x 10 _B 


0.5 


0.26 x 10" 14 


4.5 x 10" e 


1 


0.12 x 10~ 14 


1.6 x 10" 4 


2 


0.14 x 10- 1U 


1.4 x 10~ 3 


5 


0.29 x 10~ 3 


3.0 x 10~ 2 


10 


4.06 x 10 2 


8.0 



As it can be observed in the last table the result of application of MPFP in the case of 
c = 10 is less satisfactory as that of PDE tool. This is due to the fact that for larger values 
of c one should consider a bigger N. In Table 5 we show the absolute error of MPFP for 
c = 10 and N > 42. 

Table 5. Improvement in the maximum absolute error due to MPFP as the number of 

functions Uk keeps increasing 



N 


Maximum absolute error of MPFP 


42 


3.89 


44 


3.25 


46 


1.81 


48 


0.81 


50 


0.41 


52 


0.10 



Thus, one can see that for N > 42 the result obtained with the aid of MPFP is more 
accurate than that given by Matlab. We stress that in the case of using MPFP a system 
of N + 1 linear algebraic equations is solved which means solution of dozens of equations 
instead of thousands required by the finite element method implemented in the PDE tool. 

We experimented also with the shape of the domain. We considered the elliptic form as 
well as a unitary disk with a triangle shaped deformation. In the first case it is possible to 
see how the maximum absolute error increases with the excentricity e, Table 6. Here in all 
cases the area of the considered ellipses was kept constant, equal to it, while the excentricity 
was being increased. 



Table 6. Maximum absolute error for different values of the excentricity of the elliptic 
domain with the area of the domain being equal to it. The case e = corresponds to the 

unitary disc. 



TV 


e = 


e = 0.5 


e = 0.7 


e = 0.9 


e = 0.95 


e = 0.99 


30 


2.2 x 10" 14 


0.4 x 10" ia 


0.5 x 10" i3 


0.3 x 10" 12 


0.2 x 10- n 


1 x 10~ iU 



For the case of the domain with a triangular deformation (see Fig. 1), the errors were 
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Figure 1: The unitary disk with a triangular deformation. 



tested for different heights of the peak and different values of c and N with satisfactory 
results. In Table 7 we present the maximum absolute error of the approximate solution of 
the boundary value problem in dependence on the height of the triangular peak over the 
unitary circunference. 

Table 7. Maximum absolute error for N = 31, c = 1 and different heights of the peak 



Height of the peak over the unitary disk 


Maximum absolute error 


0.5 


0.92 x lO" 1 ^ 


0.7 


0.62 x 10" 11 


1.0 


0.72 x 10~ 1U 




5.2 Results obtained with numerically calculated base functions 

The use of the numerically calculated base functions which we denote by Uk poses the natural 



question about the accuracy of their calculation. Consideration of equation (31 ) gives us the 
possibility to compare u\. with the symbolically calculated exact solutions (32). In the 
following table we give the difference between ut and for c 



1. 



Table 8. Maximum absolute error of calculation of the 



aase functions for c = 1 



k 


N - u k \ 


1 


0.00000667646050 x 10" 


-4 


5 


0.00022726096338 x 10" 


-4 


11 


0.01492959925020 x 10" 


-4 


16 


0.07644765777970 x 10" 


-4 


20 


0.22545767650040 x 10" 


-4 
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Regarding the results given in the last table it is important to note that in fact the weight 
of Uk in the expansion of a solution decreases as (^p) ! or (|) ! for odd or even k respectively. 



This is due to the factor 1/n! in the definition of the Taylor coefficients (19). That is in fact 



the real accuracy of calculation of the base functions would be given by \uk — v-k \ /K\ where 
r (Ml)! f or k odd 

K = < v d\ ! ' , . It is easy to see that in the case of the results given in Table 7 

\ (I)! for k even 

one then obtains |u20 — U2o\ /10! — 6.213 x 1CP 12 which is a remarkably good agreement. 
For the numerical computation of we implemented the following procedure. Before 



integrating on each new step according to (16) along segments joining the center of the 



formal powers with points on the boundary of the domain the integrand was represented as 
a cubic spline which then was integrated using the standard Matlab routine for integration 
of splines. This procedure is simple but clearly not optimal. Nevertheless the approximate 
results presented in this work show that even such integration procedure gives satisfactory 
agreement between the exact base functions and those calculated numerically. 

The accuracy of the approximate solution obtained with the aid of in our numerical 
tests did not differ significantly from that of the solution obtained using the exactly calculated 
Uk. The order of the maximum absolute error for a given N and c coincided in both cases. 
Hence here we present results corresponding to another test problem for which we did not 
have the exactly calculated base functions. 

Consider the equation 



A + —\u(x,y) = 0. (34) 
An exact solution for this equation can be found using the fact that the change of variables 



y y 

£ = eacosf, rj = ezsinf leads to (|31|) in the new variables. Thus, e.g., the function 

y 



u{x,y) = exp(e2 cos |) is an exact solution of (34). Consequently, as a test problem we can 



consider the problem of finding a solution of (34) in $7, satisfying the boundary condition 



y X 

u(x, y) = exp(e2 cos — ), (x,y) £ d£l. (35) 

In order to construct a particular solution no in a separable form we solve numerically the 
ordinary differential equation 

d 2 e y\ . . 

v + Tr° (y)=0 - 

The obtained solution we then use for constructing the system of functions Some 
results on the accuracy of the approximate solution are given in the following table. 



Table 9. Maximum absolute error of the approximate solution of the test problem (34), 



( 35 ) considered in a unitary disk in dependence on N 
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N 


Maximum absolute error 


4 


0.079 


6 


0.021 


8 


0.005 


10 


0.001 


12 


0.0004 


14 


0.000099 


16 


0.000020 


18 


0.0000033 


20 


0.00000060 


28 


0.00000000072 


32 


0.00000000028 



6 Approximate solution of eigenvalue problems 

In this section we consider the application of MPFP to solution of eigenvalue problems for 
operators of the form —A + q(y). For simplicity we keep working with the Dirichlet boundary 
conditions and suppose that q is continuous and q(y) > in $7. Then the spectrum of the 
operator is discrete and positive. As was explained in section [4] for the equation 

(-A + q(y))u (x,y) = X 2 u {x,y) (36) 

it is easy to propose a particular solution in a separable form for any value of A. We are 
interested here in positive values, and hence a natural choice of a nonvanishing solution 



would be uo(x,y) = e lXx h(y) where h(y) is a positive solution of (28). As was observed in 
section [4] the completeness of the system of solutions {iifc}^L obtained in this case is up to 
now an open problem due to the fact that uq is complex valued and one should consider 
bicomplex pseudoanalytic formal powers for which the whole theory is still underdeveloped. 
Nevertheless we used the constructed system of exact solutions {u^^Lq for finding the eigen- 
values A 2 in the following way. Assuming that {iffcl^Lo is complete in the same sense as was 



proved in the case of the real- valued particular solution uq (subsection 3.2) we have then 



that if a nontrivial solution u of (36) exists satisfying the boundary condition u\qq = then 
u ~ X^fcLo bkUk and the coefficients bk are such that the trivial boundary condition is approx- 
imately fulfilled. This means that one can require that ^2k=o^kUk( z j) = for zj £ dQ and 
j = 0,N. This is possible iff the determinant of the matrix U = (ujk) ■ k=0 vanishes where 
u jk = Uk(zj). The determinant of U for a fixed is a function of A. Thus, the problem of 
finding eigenvalues reduces to the problem of finding zeros of the function det U(X). 

As a test problem we considered the problem of calculating the eigenvalues of the Dirichlet 
problem for the Helmholtz equation (A + A 2 ) u = 0. For every A a system of exact solutions 



{u/cl^Lo can be constructed using (32) where c should be replaced by iX. Then following 
the described scheme we looked for zeros of det U(X). As it is well known (see, e.g., j6]) the 
eigenvalues of the Dirichlet problem for the Helmholtz equation in a unitary disk are squares 
of zeros of Bessel functions J n (x). Our numerical experiments showed that a relatively small 
value of N = 21 was needed for computing the first five eigenvalues with the accuracy of 
four decimals. With A^ = 23 we obtained six first eigenvalues with the same accuracy. Thus, 
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indeed, the MPFP is clearly competitive in solving eigenvalue problems for elliptic operators. 
In applying the MPFP we detected a similar problem to that described by Alexidze in [21 
Sect. 1.13] where the method of fundamental solutions (or auxiliary sources) was applied to 
eigenvalue problems. The considered determinant shows a very fast decrement (in spite of 
this the method gives good numerical results). [2] contains references to other publications 
where different ways of using the knowledge of a system of exact solutions for numerical 
solution of eigenvalue problems were studied. In this direction further research is needed. 

7 Conclusions 

A new approach for solving boundary value and eigenvalue problems for elliptic operators in 
bounded planar domains is proposed. It is based on some classical and some new results from 
pseudoanalytic function theory which allow one to construct complete systems of solutions of 
the elliptic equations. We showed the practical applicability of the numerical method based 
on this construction, studied the rate of its convergence, accuracy and other parameters of 
its performance. 
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